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This  paper  applies  a  Bayesian  probabilistic  inferential  methodology  for  the  reconstruction  of  the  location  and  emission  rate 
from  an  actual  contaminant  source  (emission  from  the  Chalk  River  Laboratories  medical  isotope  production  facility)  using  a 
small  number  of  activity  concentration  measurements  of  a  noble  gas  (Xenon-133)  obtained  from  three  stations  that  form  part 
of  the  International  Monitoring  System  radionuclide  network.  The  sampling  of  the  resulting  posterior  distribution  of  the  source 
parameters  is  undertaken  using  a  very  efficient  Markov  chain  Monte  Carlo  technique  that  utilizes  a  multiple-try  differential 
evolution  adaptive  Metropolis  algorithm  with  an  archive  of  past  states.  It  is  shown  that  the  principal  difficulty  in  the  reconstruction 
lay  in  the  correct  specification  of  the  model  errors  (both  scale  and  structure)  for  use  in  the  Bayesian  inferential  methodology.  In  this 
context,  two  different  measurement  models  for  incorporation  of  the  model  error  of  the  predicted  concentrations  are  considered. 
The  performance  of  both  of  these  measurement  models  with  respect  to  their  accuracy  and  precision  in  the  recovery  of  the  source 
parameters  is  compared  and  contrasted. 


1.  Introduction 

In  recent  years,  there  have  been  remarkable  advances  in 
sensor  technology  for  the  environmental  monitoring  and 
surveillance  of  chemical,  biological,  or  radiological  agents 
released  into  the  atmosphere,  either  deliberately  or  acci¬ 
dentally.  Two  examples  of  operational  monitoring  sensor 
networks  are  the  deployment  of  biological  sensor  arrays 
by  the  Department  of  Homeland  Security  in  various  cities 
across  the  United  States  as  part  of  the  Bio  Watch  program 

[1]  and  the  global  network  of  radionuclide  sensors  that  form 
part  of  the  International  Monitoring  System  deployed  under 
the  auspices  of  the  Comprehensive  Nuclear-Test-Ban  Treaty 

[2] .  In  parallel  with  this  development,  there  has  been  also 
significant  theoretical  and  computational  progress  in  the 
modeling  and  simulation  of  the  transport  and  dispersion  of 
materials  released  into  the  atmosphere,  involving  ever  more 
accurate  representations  of  the  whole  complex  of  processes 
on  the  entire  spectrum  of  scales  responsible  for  the  disper¬ 
sion  of  windborne  materials  (from  the  microscale  turbulent 


motions  in  the  atmospheric  boundary  layer  that  determine 
the  short-range  dispersion  to  the  larger-scale  motions  asso¬ 
ciated  with  various  weather  systems  within  the  whole  global 
atmosphere  that  determine  the  long-range  transport). 

Lying  at  the  nexus  of  these  two  developments  is  the 
requirement  for  ever  more  sophisticated  statistical  analysis 
tools  that  can  be  used  to  fuse  the  information  embodied 
in  the  measurements  of  agent  concentration  provided  by 
the  sensors  and  in  the  predictions  of  these  concentrations 
provided  by  the  atmospheric  dispersion  model(s).  This  data- 
driven  sensor-modeling  paradigm  will  enable  the  character¬ 
ization  of  the  unknown  source  following  event  detection  by 
the  network  of  sensors  (inverse  problem),  which  in  turn  can 
lead  to  a  greatly  improved  situational  awareness  and  to  a  more 
informed  decision-making  process  for  response. 

Inverse  source  modeling  is  a  notoriously  ill-posed  prob¬ 
lem  in  the  following  sense:  (1)  the  sparse  sampling  (in 
space)  of  and  the  presence  of  multiple  sampling  times  for  the 
concentration  in  a  dispersing  cloud  of  contaminant  may  lead 
to  difficulties  of  nonuniqueness  in  the  source  reconstruction 
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(namely,  there  may  be  many  source  configurations  that  are 
consistent  with  the  limited  number  of  concentration  data); 
(2)  the  physical  processes  of  mixing  and  diffusion  of  material 
in  the  atmosphere,  with  continuing  dilution,  smooth  the 
concentration  field  leading  to  a  progressive  loss  of  high- 
frequency  information  that  could  potentially  lead  to  an 
instability  in  the  source  reconstruction;  and  (3)  actual  mea¬ 
surements  of  concentration  by  sensors  and  predictions  of 
concentrations  by  atmospheric  dispersion  model(s)  are  never 
exact  (owing  to  measurement  and  model  errors,  resp.)  which 
exacerbates  the  nonuniqueness  associated  with  incomplete 
concentration  observations  and  the  instability  associated 
with  the  intrinsic  smoothing  of  the  concentration  field  by  the 
natural  mixing  processes  in  the  atmosphere. 

To  overcome  the  problems  of  nonuniqueness  and  insta¬ 
bility  in  the  inverse  source  modeling  problem,  the  most  popu¬ 
lar  approach  is  regularization  which  emphasizes  the  selection 
of  a  “well-behaved”  source  configuration  from  the  possibly 
infinite  set  of  such  configurations  that  are  consistent  with 
the  limited  number  of  noisy  concentration  measurements. 
In  this  approach,  the  idea  is  to  find  a  source  configuration 
that  minimizes  (optimizes)  a  penalty  functional  that  involves 
two  independent  terms,  namely,  a  regularization  functional 
(source  model  norm)  that  imposes  some  constraint  on  the 
solution  and  a  misfit  functional  (also  some  form  of  norm) 
that  quantifies  the  discrepancy  between  the  measured  and 
predicted  concentrations.  For  the  application  of  regulariza¬ 
tion  to  inverse  source  modeling  problems,  see,  for  example, 
Robertson  and  Persson  [3],  Robertson  and  Langner  [4], 
Thomson  et  al.  [5],  Bocquet  [6],  and  Allen  et  al.  [7]. 

Regularization  provides  a  single  “optimal”  selection  for 
the  unknown  source  distribution  ( deterministic  solution) 
with  no  rigorous  quantification  of  the  uncertainty  in  this  par¬ 
ticular  selection.  An  alternative  to  obtaining  a  single  optimal 
source  distribution  is  to  characterize  an  ensemble  of  source 
distributions  using  the  Bayesian  statistical  paradigm  (pro¬ 
viding  a  fully  probabilistic  solution  allowing  the  uncertainty 
in  the  source  reconstruction  to  be  rigorously  assessed).  The 
Bayesian  inferential  methodology  for  source  reconstruction 
provides  fully  probabilistic  information  on  all  the  parameters 
used  to  describe  the  unknown  source  distribution.  The 
probabilistic  approach  using  a  Bayesian  inferential  scheme 
for  source  reconstruction  has  been  developed  by  Keats  et  al. 
[8],  Chow  et  al.  [9],  Senocak  et  al.  [10],  Yee  et  al.  [11],  and  Yee 
[12-15].  The  primary  limiting  factor  in  the  application  of  the 
Bayesian  inferential  methodology  for  source  reconstruction 
is  computational  and,  as  a  result,  it  has  not  been  as  widely 
used  as  the  regularization  methodology. 

Most  of  the  applications  of  the  Bayesian  inferential 
methodology  for  source  reconstruction  have  used  high- 
quality  concentration  data  from  well-designed  atmospheric 
dispersion  experiments  to  validate  the  schema.  The  objective 
of  this  paper  is  to  use  concentration  data  obtained  from  an 
operational  network  of  sensors  (more  specifically  from  a  very 
small  subset  of  the  global  network  of  radionuclide  sensors 
that  form  part  of  the  International  Monitoring  System)  to 
provide  a  real-world  test  of  source  reconstruction  based 
on  the  Bayesian  statistical  paradigm  applied  to  long-range 
atmospheric  transport  on  a  continental  or  hemispheric  scale. 


It  will  be  demonstrated  that  one  of  the  key  problems  that 
need  to  be  properly  addressed  for  the  successful  application 
of  Bayesian  inference  for  source  reconstruction  when  using 
operational  concentration  data  is  the  “correct”  specification 
of  the  expected  model  errors  arising  from  the  putative 
accuracy  of  a  long-range  forecast  (or,  alternatively,  reanalysis) 
of  meteorological  fields  and  the  concomitant  prediction 
of  material  dispersion  based  on  this  forecast  or  reanalysis 
(which  can  result  in  significant  model  errors  with  a  complex 
structure  in  the  predicted  concentrations  required  for  the 
source  inversion). 

2.  Bayesian  Probability  Theory 

In  a  remarkable  paper,  Cox  [16]  demonstrated  that  proba¬ 
bility  theory,  when  interpreted  as  logic,  is  the  only  calculus 
that  conforms  to  a  consistent  theory  of  inference.  This 
demonstration  provides  the  firm  logical  basis  for  asserting 
that  probability  calculus  is  the  unique  quantitative  theory 
of  inference.  More  specifically,  the  cornerstone  of  logical 
inference  is  embodied  by  Bayes  theorem  which  itself  is 
nothing  more  than  the  product  law  of  probability  calculus  (or 
Bayesian  probability  theory): 

p(d\l)p(d\ej)  =  p(d\l)p(d\dj).  (1) 

In  (1),  p(-)  denotes  the  probability  of  a  proposition  or 
hypothesis,  “|”  denotes  “conditional  upon”  or  “given  that”,  6 
denotes  a  set  of  parameters  of  a  model  under  consideration, 
d  denotes  the  data  used  for  the  inference,  and  I  denotes  the 
background  assumptions  or  contextual  information  available 
in  the  problem.  More  specifically,  as  applied  to  the  problem 
of  source  reconstruction,  6  can  be  identified  with  the  set 
of  parameters  used  to  characterize  the  source  distribution 
(e.g.,  location  and  emission  rate);  d  can  be  associated  with 
the  measured  concentration  data;  and  I  corresponds  to  the 
contextual  information  that  is  relevant  to  the  source  recon¬ 
struction  problem  (e.g.,  reanalyzed  or  forecast  meteorology, 
atmospheric  dispersion  model  used  to  define  the  source- 
receptor  relationship). 

The  input  quantities  for  Bayesian  inference  are  on  the 
left-hand  side  of  (1)  and  are  as  follows:  p(0  \  I)  is  the  prior 
probability  for  a  hypothesis  about  the  values  of  the  source 
parameter  vector  6  which  encodes  our  state  of  knowledge 
about  these  parameters  before  the  receipt  of  the  concentra¬ 
tion  data  d  and  p(d  \  6 , 1)  is  the  likelihood  function  and 
is  considered  to  be  a  function  of  6  for  fixed  data  d.  The 
likelihood  function  incorporates  the  information  provided 
by  the  measured  concentration  data  d  into  the  inferential 
scheme. 

The  output  quantities  provided  by  the  Bayesian  inference 
are  on  the  right-hand  side  of  (1)  and  are  as  follows:  p(d  \  I)  is 
referred  to  as  the  evidence  or  global  likelihood  and  p(6  |  d,  I) 
is  the  posterior  probability  for  the  hypothesis  about  the 
values  of  the  source  parameter  vector  6 ,  evaluated  in  light 
of  the  additional  information  provided  by  the  measured 
concentration  data  d.  The  evidence  is  given  by  the  following 
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multidimensional  integral  over  the  source  parameter  space 
[cf.  (1)]: 

p(d\i)  =  ^p(d\i)p(d\0,i)de,  (2) 

which  ensures  the  proper  normalization  of  the  posterior 
distribution  p(0  |  d,  I)  (namely,  the  condition  that  the 

integral  of  the  posterior  distribution  over  its  domain  of 
definition  is  unity).  In  the  context  of  the  determination  of  the 
plausible  values  for  the  source  parameter  vector  (parameter 
estimation  problem),  it  should  be  noted  that  the  evidence  is 
simply  a  normalization  constant  that  is  independent  of  the 
source  parameter  vector  and,  as  a  consequence,  can  therefore 
be  ignored.  However,  it  should  be  stressed  that  the  evidence 
assumes  a  central  role  in  the  problem  of  model  selection 
as  discussed  in  Yee  [15]  and,  in  this  context,  this  quantity 
cannot  be  ignored  as  in  the  parameter  estimation  problem 
considered  herein. 

In  view  of  the  fact  that  the  p(d  \  I)  is  merely  a  normaliza¬ 
tion  constant  in  the  context  of  the  source  estimation  problem, 
the  key  output  of  the  Bayesian  inference  methodology  in  this 
case  is  the  posterior  probability  p(6  |  d,  I)  which  embodies  all 
the  information  about  the  unknown  source  parameters  6.  In 
particular,  p(6  |  d,  I)  provides  the  full  solution  for  the  inverse 
source  determination  (or  source  term  estimation)  problem. 
Determining  this  quantity  requires  assigning  appropriate 
functional  forms  for  the  two  input  quantities  for  Bayesian 
inference:  (1)  the  prior  probability  p(d  \  I)  and  (2)  the  likeli¬ 
hood  function  p(d  \  6 , I). 


by  far  the  most  dominant  contribution  to  the  total  error  ej 
and  is  also  the  most  difficult  to  characterize. 

In  our  current  application,  the  model  error  arises  from 
three  primary  sources.  These  are  as  follows: 

(1)  uncertainties  in  the  representation  of  various  physical 
processes  in  the  dispersion  model; 

(2)  uncertainties  in  the  input  meteorological  fields  (initial 
and  boundary  conditions)  used  to  “drive”  the  dis¬ 
persion  model  (either  numerical  weather  prediction 
uncertainties  if  these  fields  are  obtained  as  a  forecast 
or  data  assimilation  uncertainties  for  the  state  of  the 
atmosphere  if  these  fields  are  obtained  through  a 
reanalysis); 

(3)  uncertainties  in  the  numerical  solution  of  the  model 
equations  that  characterize  the  dispersion  model 
which  comprise  both  discretization  (including  rep- 
resentativity)  errors  and  statistical  model  errors,  the 
latter  of  which  arises  from  using  necessarily  a  finite 
number  of  “marked”  fluid  particles  to  estimate  the 
mean  concentration  field  in  the  case  of  a  Lagrangian 
stochastic  model  of  dispersion. 

As  a  consequence  of  the  complexity  in  structure  of  the 
error  ej  [cf.  (3)]  for  our  current  application,  it  is  extremely 
difficult  (if  not  insuperable)  to  specify  a  priori  an  exact  value 
Oj  for  the  standard  deviation  of  ej  (J  =  1, 2, . . . ,  N).  If  the 
standard  deviation  Oj  of  was  exactly  known,  then  it  can  be 
shown  by  application  of  the  principle  of  maximum  entropy 
that  a  Gaussian  distribution  of  the  form 


2.1.  Assignment  of  the  Likelihood  Function.  To  assign  a  spe¬ 
cific  functional  form  for  the  likelihood  function,  one  needs  to 
relate  the  source  parameters  6  to  the  available  concentration 
data  d.  In  this  sense,  the  likelihood  function  defines  the 
probabilistic  model  of  how  the  concentration  data  were 
generated.  Towards  this  purpose,  if  the  concentration  data 
acquired  (by  a  sensor)  at  the  space-time  point  ( Xj,tj )  is 
represented  by  d(xj,tj ),  then  the  concentration  data  and  the 
source  parameters  6  are  assumed  to  be  related  by 

dj  =  d  (xj,  tj) 

=  C(d;xJ,tj)  +  e(xI,tJ )  (3) 

=  C,(d)  +  e,,  J  =  1,2 . N, 

where  N  is  the  total  number  of  measured  concentration  data. 
In  (3),  Cj(d)  is  the  modeled  (predicted)  mean  concentration 
at  the  /th  space-time  point.  The  predicted  concentration  is 
determined  using  an  atmospheric  dispersion  model  for  a 
source  distribution  characterized  by  the  parameter  vector  6. 

The  error  (discrepancy)  between  the  measured  concen¬ 
tration  dj  and  the  predicted  concentration  Cj(d)  is  repre¬ 
sented  symbolically  in  (3)  by  ej.  In  the  problem  considered  in 
this  paper,  there  are  two  major  contributions  to  the  error  ej, 
namely,  the  observation  or  instrument  error  in  the  measured 
concentration  dj  arising  from  the  noise  inherent  in  the  sensor 
and  the  model  error  in  the  determination  of  the  predicted 
concentration  Cj(0).  Of  these  two  errors,  the  model  error  is 


p(d|W)'nJ^“pH*2(S))’  »> 

where 


N 

x2(0)^l 


7=1 


d,  -  Cj  (0)  y 

aj  J  ’ 


(5) 


would  be  the  most  conservative  choice  for  the  direct 
probability  (or  likelihood)  of  the  concentration  data  d  = 
(d^,  . . . ,  d^f)  [17]. 

Because  Oj  is  not  known  a  priori,  it  is  useful  to  character¬ 
ize  the  uncertainty  in  the  specification  of  Oj  with  a  probability 
distribution.  Following  Yee  [15],  we  choose  this  probability 
distribution  to  be  an  inverse  gamma  distribution  with  the 
following  form: 


(p{o}  |  Sj,  a,  /3)  =  2 


r(]8) 


2  j8 


1 


aj  J  a, 


exp 

7  =  1,2,.. .,N. 


(6) 


Here,  T(x)  denotes  the  gamma  function,  a  and  /3  are  scale 
and  shape  parameters  of  the  inverse  gamma  distribution,  and 
Sj  is  the  quoted  (nominal)  estimate  for  the  true  but  unknown 
standard  deviation  Oj.  The  values  for  the  hyperparameters  a 
and  f>  are  chosen  as  a  =  7i~l  and  f  -  1  following  the  rationale 
described  in  Yee  [15]. 
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The  likelihood  function  in  (4)  and  (5)  depends  on  the 
error  standard  deviations  Oj  (/  =  1,2, ...,N)  which  are 
generally  unknown.  To  remove  these  unwanted  parameters 
(nuisance  parameters),  we  can  multiply  the  likelihood  given 
in  (4)  by  the  (assigned)  probability  distribution  for  each  of 
the  error  standard  deviations  embodied  in  (6)  and  integrate 
the  result  with  respect  to  the  unwanted  parameters  (error 
standard  deviations)  to  give  an  integrated  likelihood  function 
with  the  following  form  [15]: 

p(d  |  0,  s,a,/3) 


p  (d  |  0, 1)  Y[<P  (o>  I  Sj,  a,  P)  da 

7=1 


A«^r(j8+i/2) 

7=1  '/2nsJT(p) 


oc  + 


( dj  Cj  ( 0 )) 
2? 


-P- 1/2 


(7) 


In  (7),  s  =  (s, ,  s2, . . . ,  sN)  is  the  vector  of  estimated  (quoted) 
standard  deviations  for  the  error  er{J  =  1,2  and  da  = 

dcr  ^dcr 2  •  •  •  d(Tjsj. 

To  be  more  specific  with  respect  to  the  source  parameter 
vector  0,  a  source  distribution  S  with  the  following  form  is 
considered  in  this  paper: 


s  (x)  =  QsS  (x  -  xs) ,  (8) 


where  d(-)  is  the  Dirac  delta  function.  Equation  (8)  describes 
a  (continuous)  point  source  located  at  the  vector  position 
xs  that  is  emitting  contaminant  at  a  constant  emission  rate 
Qs.  The  parameters  describing  this  source  distribution  can 
be  assembled  into  a  source  parameter  vector  given  by  0  s 

(xs,Qs). 


2.2.  Assignment  of  the  Prior.  To  assign  the  prior  probability 
for  the  source  parameters  0,  it  is  necessary  to  state  explicitly 
what  is  known  about  these  parameters.  Firstly,  it  is  assumed 
that  the  source  parameters  are  logically  independent  with  the 
result  that  the  prior  distribution  p(6  \  I)  factorizes  as  follows: 

p(0\I)  =  p(xs\l)p(Qs\l).  (9) 

Secondly,  the  source  location  and  emission  rate  are  known 
a  priori  to  be  bounded.  In  particular,  it  is  assumed  that  the 
location  xs  of  the  source  is  contained  in  some  spatial  region 
Q)  c  R 3.  Furthermore,  the  emission  rate  Qs  is  assumed 
to  be  bounded  by  Qmin  <  Qs  <  Qmax  where  Qmin  and 
Qmax  are  the  lower  and  upper  bounds  for  the  emission  rate, 
respectively.  Finally,  if  nothing  else  is  known  about  these 
parameters  except  for  the  bounds,  then  application  of  the 
principle  of  maximum  entropy  to  our  state  of  knowledge 
concerning  the  source  parameters  results  in  the  assignment 
of  a  uniform  prior  distribution  for  these  parameters,  so 

P  (0  I  I)  oc  Ds  (xs)  »(Qmin,Qmax)  (Qs)  >  (10) 

where  QA(x)  denotes  the  indicator  function  for  set  A  defined 
as  0A(x)  =  1  if  x  e  A  and  QA(x)  =  0  if  x  i  A. 


2.3.  Posterior  Distribution:  The  Output.  The  key  output  of  the 
Bayesian  inference  is  the  posterior  distribution  p(6  |  d,  I) 
which  can  be  obtained  by  combining  the  assignments  for  the 
likelihood  function  and  the  prior  distribution  given  by  (7) 
and  (10),  respectively,  to  give  [cf.  (1)] 

p(0  |  d,  s,  a,  p,  I) 

a  n3  (xs)  n{Qmln,Qmax)  (Qs) 
x  AapT(p+ 1/2)  (11) 

j=i  ^/2nsjT(p ) 

x  [a  +  (dj  -  C;(0))2/(2s^)]  P  1/2. 

Note  that  the  quantities  s,  a,  and  f>  have  been  added  explicitly 
to  the  posterior  probability  of  the  source  parameters  0  in  (11) 
to  indicate  that  these  quantities  are  known  (namely,  they  are 
provided  a  priori  by  the  user,  in  addition  to  the  measured 
concentration  data  d). 


3.  Computational  Aspects 

There  are  two  computational  problems  that  need  to  be 
addressed  before  the  Bayesian  inferential  methodology  for 
source  reconstruction  can  be  applied  to  practical  real-world 
applications;  namely,  (1)  that  Bayesian  inversion  of  concen¬ 
tration  data  requires  a  fast  and  efficient  technique  for  the 
determination  of  the  source-receptor  relationship  (namely, 
for  the  fast  computation  of  Cj(0)  for  a  given  hypothesis 
about  the  source  parameters  0)  and  (2)  methodology  for 
the  efficient  sampling  of  the  posterior  distribution  p(6  \ 
d,  s,  a,  j3, 1). 

To  address  the  first  problem,  Keats  et  al.  [8]  and  Yee 
et  al.  [11]  demonstrated  that  the  use  of  a  receptor-oriented 
representation  for  the  source-receptor  relationship  (rather 
than  the  more  usual  source-oriented  representation)  enabled 
the  efficient  computation  of  the  predicted  concentration  at 
a  fixed  receptor  location  associated  with  the  exploration  of 
a  large  number  of  source  parameter  hypotheses  required 
in  the  simulation-based  Bayesian  inference  methodology. 
For  the  application  reported  herein,  a  backward  Lagrangian 
stochastic  (LS)  model  for  long-range  transport  was  used 
to  determine  the  adjoint  concentration  field  C*  over  the 
northern  hemisphere. 

The  backward  LS  model  employed  here  is  an  opera¬ 
tional  model  used  by  the  Canadian  Meteorological  Centre 
to  support  both  Canadian  Treaty  monitoring  and  the  var¬ 
ious  mandates  of  the  Provisional  Technical  Secretariat  of 
the  Comprehensive  Nuclear-Test-Ban  Treaty  Organization, 
including  event  analysis  by  member  states.  The  backward 
LS  model  for  the  determination  of  C*  was  “driven”  by 
reanalyzed  meteorological  fields  that  were  obtained  at  a 
relatively  low  temporal  and  spatial  resolution,  namely,  at  a 
temporal  resolution  of  6  h  (rate  of  the  data  assimilation)  with 
a  core  spatial  resolution  of  0.5°  on  a  geographical  latitude  and 
longitude  coordinate  system.  The  backward  LS  model  was  run 
retrospectively  using  these  reanalyzed  meteorological  fields, 
providing  C*  fields  with  a  temporal  resolution  of  3  hours  over 
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a  period  of  14  days  prior  to  the  commencement  of  the  sam¬ 
pling  for  a  particular  activity  concentration  measurement. 

It  should  be  noted  that,  with  the  source  distribution 
S  given  by  (8),  the  predicted  concentration  at  a  receptor 
location  Xj  and  time  tj  can  be  easily  determined  from  the 
following  relationship  once  the  adjoint  concentration  field  C* 
has  been  computed  for  this  location  and  time: 

Cj(d)  =  c(d-,xj,tj) 

=  [  dt'  [  dxC*  (x,t'  |  Xr,t,)s(x)  (12) 
J-oo  Js 5  v  /  V  / 

®  (C  I  S)  (xj,  tj ) , 

where  C*(x\  t'  |  Xj,  tj )  is  the  adjoint  concentration  at  space- 
time  point  (x,tr)  associated  with  the  sensor  concentration 
datum  measured  at  location  Xj  and  time  tj.  In  (12),  is 
a  volume  in  space  that  contains  the  source.  The  reader  is 
referred  to  Marchuk  [18]  and  Le  Dimet  and  Talagrand  [19] 
for  a  more  detailed  explanation  of  the  adjoint  formulation  in 
a  more  general  context. 

Note  that  the  predicted  mean  concentration  Cj(6)  “seen” 
by  a  sensor  for  a  given  hypothesis  about  the  source  parame¬ 
ters  6  can  be  rapidly  computed  by  simply  evaluating  the  inner 
(or  scalar)  product  (C*  |  S)  of  the  adjoint  concentration 
C*  and  the  source  distribution  S  corresponding  to  the  given 
hypothesis.  In  other  words,  the  predicted  concentration  Cj  is 
obtained  from  a  mathematical  model  (backward  LS  model  for 
dispersion)  by  evaluation  of  the  bounded  linear  functionals 
Cj  =  (Cj  |  S)  for  /  =  1,2, ...,N  where  C*j  denotes  the 
adjoint  concentration  field  obtained  at  the  sensor  space- time 
point  ( Xj,tj )  and  where  S  corresponds  to  a  given  hypothesis 
about  the  source.  In  applied  mathematics,  the  elements  Cj 
are  referred  to  usually  as  representers.  More  specifically,  if  we 
substitute  (8)  into  (12),  the  model  (predicted)  concentration 
Cj(0)  “seen”  by  the  sensor  at  location  Xj  and  time  tj  is  given 
explicitly  by 

Cj  (6)  =  Qs  [  C*  (xs,  t'  I  Xj,  tj )  dt'.  (13) 

J-OO 

To  address  the  second  computational  problem  mentioned 
above,  we  apply  a  Markov  chain  Monte  Carlo  (MCMC)  algo¬ 
rithm  for  the  efficient  sampling  of  the  posterior  distribution 
of  the  source  parameters  (see  Gilks  et  al.  [20],  Gelman  et  al. 
[21],  Savchuk  and  Tsokos  [22],  and  Yuen  [23]).  In  general, 
MCMC  algorithms  generate  a  Markov  chain  whose  station¬ 
ary  distribution  coincides  exactly  with  the  target  probability 
distribution  that  we  are  trying  to  sample  from  (which  in 
our  case  is  the  posterior  distribution  p(6  |  d,  s,  a,  /3 , 1)  given 
in  (11)).  The  Metropolis-Hastings  (M-H)  algorithm  [20,  21] 
forms  the  underlying  basis  for  MCMC  sampling  and  it  is 
perhaps  not  too  surprising  that  the  M-H  algorithm  has 
become  almost  synonymous  with  MCMC  sampling.  Indeed, 
most  of  the  algorithms  for  MCMC  sampling  reported  in  the 
literature  [20-23]  can  be  interpreted  as  either  special  cases  or 
extensions  of  the  basic  M-H  algorithm. 

For  the  current  study,  we  use  a  multiple-try  differ¬ 
ential  evolution  adaptive  Metropolis  algorithm  with  sam¬ 
pling  from  an  archive  of  past  states  which  is  referred  to 


as  MT-DREAM(ZS).  The  details  of  this  MCMC  sampling 
algorithm  are  described  by  Laloy  and  Vrugt  [24],  but,  for 
completeness,  we  will  briefly  summarize  the  main  compo¬ 
nents  of  this  algorithm.  In  particular,  only  the  relevant  details 
of  the  algorithm  that  are  required  for  the  interpretation 
of  the  results  in  this  paper  are  emphasized.  Firstly,  MT- 
DREAM(ZS)  samples  from  an  archive  of  past  states  to  generate 
the  candidate  points  (proposals)  for  each  of  the  individual 
Markov  chains  that  are  used  to  explore  the  target  posterior 
distribution.  Secondly,  as  already  alluded  to  here,  the  algo¬ 
rithm  utilizes  multiple  (different)  Markov  chains  that  are 
run  simultaneously  in  parallel.  These  multiple  chains  employ 
a  self-adaptive  randomized  subspace  sampling  of  difference 
vectors  from  the  archive  of  past  states  to  generate  new 
candidate  points  in  these  chains.  As  part  of  the  randomized 
subspace  sampling  strategy,  each  element  of  the  candidate 
points  for  the  parallel  proposals  is  updated  (accepted)  in 
accordance  with  a  binomial  scheme  (Bernoulli  trial)  with 
a  crossover  probability  ps;  otherwise,  the  proposed  element 
retains  its  previous  (old)  value.  This  multiple -chain  approach 
automatically  adjusts  or  adapts  the  scale  and  orientation  of 
the  proposal  function.  Thirdly,  a  snooker  step  update  of  the 
state  with  an  adaptive  step  size  is  included  with  a  fixed  (albeit 
small)  probability  in  order  to  improve  the  mixing  efficiency 
of  the  algorithm  for  exploration  of  the  hypothesis  space.  The 
updated  states  from  the  various  parallel  Markov  chains  are 
periodically  stored  in  the  archive  of  past  states  after  every  K 
(K  >  1)  updates. 

Fourthly,  to  further  improve  the  efficiency  of  the  sam¬ 
pling,  MT-DREAM(ZS)  incorporates  a  multiple-try  Metropo¬ 
lis  (MTM)  approach  proposed  initially  by  Liu  et  al.  [25].  The 
basic  idea  underpinning  the  MTM  approach  is  as  follows: 
longer  range  candidate  moves  are  rarely  accepted,  but  if 
multiple  points  are  proposed  for  these  longer  range  moves 
then  the  acceptance  probability  will  be  increased.  The  MTM 
algorithm  is  applied  individually  to  each  of  the  different 
Markov  chains  used  in  the  MT-DREAM(ZS),  involving  gen¬ 
erating  K  draws  using  the  randomized  subspace  sampling 
procedure  for  each  chain,  choosing  one  of  these  draws 
(proposals)  as  the  reference  point,  and  generating  a  new  set 
of  (K  -  1)  draws  with  respect  to  this  reference  point  using 
the  randomized  subspace  sampling  strategy.  The  acceptance 
rule  is  simply  the  Metropolis-Hasting  acceptance  probability 
applied  to  the  sequence  of  proposals  that  comprise  the  MTM 
schema. 

Finally,  to  determine  if  the  multiple  Markov  chains  used 
in  MT-DREAM(ZS)  have  achieved  stationarity  (namely,  have 
converged  to  the  stationary  distribution  associated  with  the 
chains),  the  Gelman  and  Rubin  [26]  convergence  diagnostic 
R  is  computed  for  each  dimension  of  each  chain  using  the 
last  50%  of  the  samples  of  the  chain.  This  simple  convergence 
statistic  determines  whether  a  chain  has  achieved  stationarity 
by  comparing  the  variance  within  each  chain  and  the  variance 
between  chains  (interchain  variance). 

4.  Example:  A  Real-World  Application 

The  International  Monitoring  System  (IMS)  consists 
of  a  comprehensive  network  of  seismic,  hydroacoustic, 
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Figure  1:  Locations  of  the  three  sampling  stations  from  the  Interna¬ 
tional  Monitoring  System  radionuclide  network  used  for  the  source 
reconstruction.  The  release  location  of  the  Xe-133  tracer  source  was 
at  Chalk  River  Laboratories. 


infrasound,  and  radionuclide  sensors  as  part  of  the 
verification  regime  of  the  Comprehensive  Nuclear-Test- 
Ban  Treaty  (CTBT)  which  bans  nuclear  explosions.  A 
subset  of  the  IMS  is  the  subnetwork  of  radionuclide  gamma 
detectors/particle  filters  for  the  measurement  of  the  activity 
concentration  for  various  radionuclides  (e.g.,  particulate 
or  aerosol  species  such  as  Cs-137  and  1-131  and/or  noble 
gases  such  as  Xe-133).  The  IMS  radionuclide  network 
will  (eventually)  have  80  monitoring  stations  worldwide 
for  the  measurement  of  the  activity  concentration  for 
particulate/aerosol  radioactive  species,  of  which  at  least 
40  stations  would  also  have  the  capability  to  measure  the 
activity  concentration  of  noble  gases  [2] .  The  stations  provide 
12  or  24  h  averaged  activity  concentrations  of  the  various 
radionuclides  depending  on  the  technology  used. 

We  applied  the  proposed  source  reconstruction  method¬ 
ology  to  some  (albeit  limited)  measurements  of  activity 
concentrations  by  a  very  small  subset  of  the  IMS  radionuclide 
network.  The  release  event  consisted  of  the  stack  emissions 
from  Chalk  River  Laboratories  (CRL)  which  houses  an 
international  production  facility  for  medical  radioisotopes. 
Chalk  River  Laboratories,  which  at  a  latitude  of  46.15°N  and  at 
a  longitude  of-77.37°E,  is  located  about  180  km  northwest  of 
the  city  of  Ottawa,  Ontario.  A  characterization  of  the  weekly 
stack  emissions  of  Xe-133  from  the  CRL  medical  isotope 
production  facility  over  a  5 -year  period  yielded  a  median 
daily  emission  of  about  24  TBq  (or,  equivalently,  an  emission 
rate  of  about  1.0  x  1012  Bqh-1). 

Xe-133  activity  concentrations  were  detected  and  mea¬ 
sured  by  three  sampling  sites  in  North  America.  These  three 
sampling  sites  are  shown  in  Figure  1.  The  three  sites  displayed 
are  part  of  the  noble  gas  monitoring  network  of  the  CTBT 
verification  regime.  More  specifically,  the  three  stations  are  as 
follows:  CAX17  (St.  Johns,  Newfoundland,  at  latitude  47.59° N 
and  longitude  -52.74°E);  USX75  (Charlottesville,  Virginia,  at 
latitude  38.0°N  and  longitude  -78.4°E);  and  USX74  (Ashland, 
Kansas,  at  latitude  31.17°N  and  longitude  -99.77°E). 


At  each  monitoring  site,  one  of  two  different  moni¬ 
toring  technologies  is  employed  to  measure  radioxenon. 
The  St.  Johns  site  has  a  Systeme  de  Prelevements  et 
dAnalyse  en  Ligne  dAir  pour  Quantifier  le  Xenon  (SPALAX) 
high- resolution  gamma  system  operating  on  a  24  h  sample 
collection  period,  while  the  remaining  sites  have  a  Swedish 
Automatic  Unit  for  Noble  Gas  Acquisition  (SAUNA)  beta- 
gamma  coincidence  system  with  a  12  h  sample  collection 
period.  Both  systems  employ  activated  charcoal  to  remove 
xenon  from  the  air  for  radioxenon  analysis.  After  the  mea¬ 
surement  process  is  complete,  the  xenon  sample  can  be  stored 
in  an  archive  bottle  for  optional  remeasurement  either  on¬ 
site  or  in  an  off- site  laboratory.  The  stable  xenon  volume 
collected  varies  approximately  from  2  mL  to  8  mL  depending 
on  the  technology  used,  but  both  technologies  have  a  roughly 
equivalent  Xe-133  sensitivity  of  approximately  0.2  mBq  m-3. 

The  monitoring  data  for  activity  concentrations  of  Xe- 
133  were  obtained  from  a  single  month  (December  2011).  The 
three  monitoring  stations  used  for  this  study  were  known 
to  be  operating  normally  during  this  period  of  time.  For 
the  purposes  of  source  reconstruction,  we  used  36  of  these 
concentration  samples  extracted  from  the  three  sampling 
sites.  The  concentration  time  series  data  at  each  station  were 
blocked  averaged  over  a  temporal  duration  of  3  days  to  give 
8  concentration  data  dj  that  were  used  as  the  input  for 
the  source  reconstruction  methodology  (namely,  N  -  8 
for  this  example).  It  should  be  noted  that  predictions  for 
the  temporally  averaged  concentration  data  were  subject  to 
smaller  statistical  sampling  errors  (for  a  fixed  number  of 
marked  fluid  particles  used  in  the  backward  LS  model), 
which  is  the  primary  reason  for  aggregating  the  12  or 
24  h  measurements  in  this  manner.  As  mentioned  earlier, 
reanalyzed  meteorological  fields  were  used  to  “drive”  an 
operational  backward  LS  model  to  determine  C*  which  were 
utilized  subsequently  for  the  fast  calculation  of  the  predicted 
concentration  Cj(0)  for  an  arbitrary  hypothesis  6  about  the 
source. 

In  addition  to  the  measured  activity  concentration  data 
dj  and  the  predicted  activity  concentration  data  Cj(0 ),  it 
is  also  necessary  to  provide  an  estimate  for  the  noise  error 
variance  s2  [cf.  (7)].  The  estimate  for  the  noise  error  variance 
includes  two  major  contributions;  namely,  (1)  an  estimate  for 
the  sensor  error  sampling  variance  s2  j  in  the  measurement 
of  dj  and  (2)  an  estimate  for  the  model  error  variance 
s2m  J  in  the  prediction  of  Cj(6).  These  estimates  for  the  two 
contributions  to  the  noise  error  variance  were  added  in 
quadrature  to  give  s2  =  s2eJ  +  s2m  J.  Very  good  estimates 
for  the  sensor  error  standard  deviation  (or  square  root  of 
the  variance)  were  provided  for  the  expected  precision  in 
the  measurements  of  the  activity  concentration  at  each  of 
the  three  sampling  stations.  In  contrast,  no  estimates  for  the 
precision  in  the  predicted  concentrations  were  available.  In 
light  of  this,  the  model  error  standard  deviation  was  assumed 
(rather  arbitrarily)  to  be  50%  of  the  predicted  concentration 
Cj(0). 

For  our  example,  the  emission  source  is  assumed  to  be  at 
or  near  ground  level  so  that  the  height  of  the  source  is  not  of 
any  interest  in  the  reconstruction.  As  a  result,  the  unknown 
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Figure  2:  Univariate  (diagonal)  and  bivariate  (off- diagonal)  marginal  posterior  distributions  of  the  source  parameters;  namely,  longitudinal 
position  xgy  latitudinal  position  ys,  and  emission  rate  Qs.  The  true  parameter  values  are  shown  by  a  solid  square  or  a  solid  vertical  line  and 
the  best  estimates  of  the  parameter  values  are  represented  as  a  solid  circle  or  a  dashed  vertical  line. 


location  parameters  for  the  source  are  its  longitudinal  (xs) 
and  latitudinal  (ys)  positions,  so  6  =  (xs>ys>Qs).  The 

proposed  stochastic  sampling  algorithm  constructs  an  initial 
population  for  the  archive  of  past  states  by  sampling  from  the 
prior  distribution  p(6  \  I)  (see  (10)).  The  hyperparameters 
defining  the  prior  distribution  p(0  \  I)  are  chosen  as  follows: 
Qmin  =  1.0  X  1015f<Bqh-1,  Qmax  =  1.0  x  lO^Bqh-1 
which  prescribes  the  lower  and  upper  bounds  for  the  source 
emission  rate  and  2d  =  (-125°E, -45°E)  x  (25°N,  75°N) 
which  defines  the  prior  bounds  for  the  location  (xs>  ys)  of 
the  source.  It  is  noted  that  2d  encloses  the  North  American 
continent. 

The  samples  of  6  drawn  from  the  posterior  distribution 
p(0  |  d,  s,  a,  /3, 1)  (cf.  (11))  using  the  MT-DREAM(ZS)  sam¬ 
pling  algorithm  were  used  to  determine  the  characteristics 
(location  and  emission  rate)  of  the  source.  Estimates  of  these 
source  parameters  were  obtained  directly  from  the  multiple 
chains  of  samples  generated  by  the  MCMC  algorithm  after 
convergence  of  the  chains  as  determined  by  the  Gelman- 
Rubin  R  statistic.  Figure  2  exhibits  the  univariate  (diagonal) 
and  bivariate  (off-diagonal)  marginal  posterior  distributions 
for  the  source  parameters.  For  the  univariate  marginal 


posterior  distribution  of  a  parameter,  the  solid  vertical  line 
delineates  the  true  value  of  the  parameter,  and  the  dashed 
vertical  line  corresponds  to  the  best  estimate  of  the  parameter 
obtained  as  the  posterior  mean.  Similarly,  for  the  bivariate 
marginal  posterior  distribution  of  various  combinations  of 
parameters,  the  solid  square  marks  the  position  of  the  true 
source  parameter  values  and  the  solid  circle  exhibits  the  best 
estimate  of  the  true  source  parameter  values  obtained  as 
the  posterior  mean.  It  should  be  noted  that  the  axes  limits 
are  chosen  to  display  the  regions  of  the  highest  probability 
in  the  marginal  posterior  probability  distributions  for  the 
parameters  and,  in  certain  cases,  these  regions  do  not  contain 
the  true  parameter  values  (with  the  result  that  some  of  the 
panels  in  Figure  2  do  not  contain  either  the  solid  square  or 
solid  vertical  line  representing  the  true  parameter  values). 

The  posterior  mean,  posterior  standard  deviation,  and 
lower  and  upper  bounds  for  the  95%  highest  posterior 
distribution  (HPD)  interval  (the  p%  HPD  (or  credible) 
interval  is  that  interval  that  contains  the  parameter  with  a  p% 
(posterior)  probability,  with  the  lower  and  upper  bounds  of 
the  interval  specified  such  that  the  probability  density  within 
the  interval  is  everywhere  larger  than  that  outside  it)  for  the 
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Table  1:  The  posterior  mean,  posterior  standard  deviation,  and  lower  and  upper  bounds  of  the  95%  HPD  interval  of  the  parameters  (°E), 
ys  (°N),  and  Qs  (pBq  IT1)  calculated  from  samples  of  0  drawn  from  the  posterior  distribution  p(0 |d,  s,  a,  1). 


Parameter 

Mean 

Standard  deviation 

95%  HPD 

Actual 

(°E) 

-72.71 

1.22 

(-76.71,  -71.56) 

-77.37 

CN) 

42.51 

0.32 

(41.86,  43.35) 

46.15 

Qs  (f'Bqh-1) 

6.68  x  1017 

8.55  x  1016 

(4.98,  8.33)  xlO17 

1.0  x  1018 

Figure  3:  Two-dimensional  marginal  posterior  distribution  of  the 
source  location  georeferenced  on  a  Google  Earth  image. 


reconstructed  source  parameters  are  summarized  in  Table  1. 
It  is  noted  that  the  accuracy  in  the  reconstruction  of  the 
source  parameters  is  fairly  good  for  both  the  location  (con¬ 
sidering  the  fact  that  the  reconstruction  was  undertaken  on 
a  continental  scale)  and  the  emission  rate.  More  specifically, 
the  distance  between  the  true  source  location  and  the  best 
estimate  of  the  source  location  (obtained  as  the  posterior 
mean)  is  only  about  572  km  (see  Figure  3)  and  the  best 
estimate  of  the  emission  rate  (obtained  again  as  the  posterior 
mean)  is  within  33%  of  the  true  emission  rate. 

Even  so,  a  perusal  of  Table  1  indicates  that  the  precision 
in  the  estimates  of  both  the  source  location  and  the  emission 
rate  is  poorly  determined.  Indeed,  it  should  be  noted  that 
the  95%  HPD  intervals  for  longitudinal/latitudinal  positions 
and  the  emission  rate  of  the  reconstructed  source  do  not 
contain  the  true  source  parameters.  This  defect  in  the  source 
reconstruction  can  be  attributed  to  the  difficulties  in  provid¬ 
ing  good  estimates  for  the  model  errors  in  the  prediction 
of  Cj(0).  It  is  evident  that  the  inability  to  provide  good 
estimates  for  the  model  errors  (which  in  principle  can  be  quite 
significant  in  our  example)  can  lead  to  a  loss  of  power  in  the 
source  reconstruction  and  may  mask  important  features  of 
the  measured  activity  concentration  data. 

Given  the  complexity  in  the  heteroskedastic  variance  of 
the  model  error  in  our  application,  the  specification  of  the 
prior  information  on  the  structure  of  the  model  error  is 
extremely  difficult.  In  light  of  these  seemingly  insuperable 
difficulties,  we  consider  an  alternative  measurement  model 


to  that  introduced  earlier  in  (3).  To  this  purpose,  let  us  now 
focus  on  the  following  alternative  measurement  model: 

dj  -  mjCj  ( 6 )  +  n7,  /  =  1, 2, . . . ,  N,  (14) 

where  m7  are  unknown  multipliers  (scale  factors)  that  are 
applied  to  the  predicted  concentration  Cj(0)  in  order  to  com¬ 
pensate  for  the  model  uncertainty.  It  is  important  to  note  that, 
in  (14),  the  rij  term  represents  only  the  measurement  error 
in  the  activity  concentration  dj.  The  model  errors  incurred 
in  the  predicted  concentration  Cj(0)  are  compensated  by  the 
introduction  of  the  multiplier  mj. 

For  the  alternative  measurement  model,  the  multipliers 
mj  (/  =  1, 2, . . . ,  N)  are  unknown  parameters  and  need 
to  be  estimated  in  addition  to  the  source  parameters.  Let 
us  denote  the  source  parameters  in  this  alternative  model 
by  6s  =  (xs,  Qs)  and  by  6m  =  (mv  m2, . . . ,  mN)  all  other 
relevant  parameters  (multipliers  in  our  example).  These  latter 
parameters  are  usually  referred  to  as  nuisance  parameters. 
Both  sets  of  parameters  define  the  parameter  vector  as 
6  =  (0s,  0m).  The  likelihood  function  for  the  alternative 
measurement  model  is  still  given  by  (7)  with  the  implicit 
understanding  that  the  estimated  noise  variances  s2  now 
only  include  the  measurement  error  contribution  (namely, 
s2  =  s2eJ).  As  already  noted  above,  this  contribution  to  the 
uncertainty  is  well  characterized  in  our  current  application 
(implying  that  the  prior  uncertainty  in  the  measurement 
errors  can  be  specified  correctly).  For  the  alternative  mea¬ 
surement  model,  one  requires  also  the  specification  of  the 
prior  distributions  for  the  multipliers  (nuisance  parameters). 
Towards  this  objective,  uniform  priors  defined  over  the 
range  from  mmin  to  ramax  will  be  adopted  as  priors  for  the 
multipliers.  In  view  of  this,  the  prior  distribution  for  the 
alternative  measurement  model  replaces  (10)  by  the  following 
assignment: 

N 

P  (0  I  I)  OC  »2>  (xs)  “(Qmin,Qmax)  (Qs)  X  n  (m/)  • 

7=1 

(15) 

Using  the  new  measurement  model  and  the  modified 
likelihood  function  and  prior  distribution,  we  applied  the 
MT-DREAM(ZS)  algorithm  to  sample  from  the  modified  pos¬ 
terior  distribution  for  0.  The  hyperparameters  that  define  the 
prior  distribution  for  the  source  parameters  0s  were  exactly 
as  described  above.  The  prior  bounds  for  the  multipliers  m7 
were  mmin  =  0.1  and  mmax  =  10.0  for  /  =  1,2,  ...,8 
(recall  that  N  =  8  in  our  example).  The  one-dimensional 
and  two-dimensional  marginal  posterior  distributions  of  the 
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Table  2:  The  posterior  mean,  posterior  standard  deviation,  and  lower  and  upper  bounds  of  the  95%  HPD  interval  of  the  parameters  (°E), 
ys  (°N),  and  Qs  (^BqfT1)  calculated  from  samples  of  0  drawn  from  the  posterior  distribution  p(0 |d,  s,  a,  /3, 1).  The  results  are  obtained  using 
the  alternative  measurement  model  which  employs  multipliers  to  compensate  for  the  unknown  model  error  structure. 
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Figure  4:  Univariate  (diagonal)  and  bivariate  (off-diagonal)  marginal  posterior  distributions  of  the  source  parameters;  namely,  longitudinal 
position  xgy  latitudinal  position  ys,  and  emission  rate  Qs.  The  true  parameter  values  are  shown  by  a  solid  square  or  a  solid  vertical  line  and 
the  best  estimates  of  the  parameter  values  are  represented  as  a  solid  circle  or  a  dashed  vertical  line.  The  reconstruction  was  obtained  using 
the  alternative  measurement  model  which  employs  multipliers  to  compensate  for  the  unknown  model  error  structure. 


source  parameters  0s  are  displayed  in  Figure  4.  The  true 
source  parameters  represented  by  either  the  solid  square  or 
solid  vertical  line  should  be  compared  with  the  best  estimates 
of  these  parameters  (posterior  mean)  marked  by  either  a 
solid  circle  or  dashed  vertical  line.  Table  2  summarizes  the 
posterior  mean,  posterior  standard  deviation,  and  lower  and 
upper  bounds  for  the  95%  HPD  interval  for  the  source 
parameters. 

It  is  clear  from  an  examination  of  Figure  4  and  of 
Table  2  that  all  the  source  parameters  have  been  recovered 


with  very  good  accuracy.  More  specifically,  the  distance 
between  the  actual  source  location  and  the  best  estimate 
(posterior  mean)  of  the  source  location  is  about  44  km  (see 
Figure  5),  and  the  recovered  emission  rate  is  within  20%  of 
the  actual  emission  rate.  Note  that  the  accuracy  in  the  inferred 
source  location  is  roughly  a  factor  of  ten  better  than  that 
obtained  using  the  standard  measurement  model  (which  does 
not  use  multipliers  to  try  to  compensate  for  the  unknown 
model  errors).  More  importantly,  the  precision  of  the  source 
parameter  estimates  (namely,  the  95%  credible  intervals)  for 
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Figure  5:  Two-dimensional  marginal  posterior  distribution  of 
the  source  location  georeferenced  on  a  Google  Earth  image.  The 
marginal  distribution  is  for  the  alternative  measurement  model 
which  employs  multipliers  to  compensate  for  the  unknown  model 
error  structure. 


the  alternative  measurement  model  which  utilizes  multipliers 
to  compensate  for  model  errors  contains  the  actual  (true) 
values  for  the  source  parameters,  in  stark  contrast  to  the 
reconstruction  using  the  standard  measurement  model  (cf. 
the  inferred  values  for  the  source  parameters  in  Table  1  with 
those  in  Table  2). 

In  the  practical  real-world  example  of  source  recon¬ 
struction  considered  herein  using  activity  concentration  data 
obtained  from  the  IMS  radionuclide  network,  it  is  extremely 
difficult  to  correctly  specify  a  priori  both  the  magnitude  (or 
scale)  of  the  model  error  and  the  model  error  structure.  In 
light  of  this  difficulty,  it  is  apparent  from  this  example  that 
incorporating  multipliers  (scale  factors)  with  the  predicted 
concentrations  to  attempt  to  compensate  for  the  model 
error  can  improve  significantly  the  quality  of  the  source 
reconstruction.  As  an  indication  of  the  complexity  in  the 
model  error  structure  for  the  current  problem,  Figure  6 
exhibits  the  one-dimensional  marginal  posterior  distribu¬ 
tions  for  the  various  multipliers  mj  (/  =  1,2,  ...,8). 

An  examination  of  Figure  6  provides  an  indication  of  the 
complexity  in  the  heteroskedastic  model  error  structure, 
which  makes  it  difficult  to  provide  a  priori  estimates  for  the 
model  error  for  use  in  the  source  reconstruction.  Note  that 
the  distributions  for  the  multipliers  associated  with  some 
of  the  predicted  concentrations  are  quite  broad,  indicating 
that  the  model  errors  for  these  predicted  concentrations  are 
not  well  determined.  In  other  cases,  the  distributions  for 
the  multipliers  (e.g.,  m2  and  m3)  are  quite  narrow,  implying 
that  the  data  contain  reasonable  information  to  constrain 
the  model  errors  associated  with  these  multipliers  fairly 
well. 

The  true  values  for  the  multipliers  are  not  known. 
However,  we  can  obtain  an  independent  estimate  for  the 
values  of  the  multipliers  in  the  current  example  because 
the  actual  source  parameters  are  known.  We  can  use  the 


actual  source  parameters  6*  to  predict  the  activity  con¬ 
centration  Cj(0*)  that  would  be  expected  at  the  various 
observation  stations.  With  this  information,  we  can  provide 
an  independent  point  estimate  for  the  multipliers  as  rhj  - 
dj/Cj(d *),  /  =  1,2, ...,8.  These  “actual”  values  for  the 

multipliers  are  shown  as  the  solid  vertical  lines  in  Figure  6, 
which  should  be  compared  with  best  estimates  of  the  mul¬ 
tipliers  obtained  as  the  mean  of  the  posterior  distributions 
for  mj.  A  visual  inspection  shows  that  the  best  estimates 
of  the  multipliers  are  consistent  with  the  “actual”  values  for 
the  multipliers  obtained  from  the  indicated  point  estimates. 
Furthermore,  the  width  of  the  posterior  distribution  of  the 
various  multipliers  is  seen  to  enclose  the  “actual”  values  for 
the  multipliers.  Finally,  even  though  some  of  the  multipliers 
are  highly  uncertain,  inclusion  of  multipliers  to  compensate 
for  the  model  errors  does  lead  to  significantly  improved 
estimates  for  the  source  parameters  (both  in  accuracy  and  in 
precision). 

5.  Conclusions 

In  this  study,  a  Bayesian  inferential  methodology  has  been 
applied  to  the  problem  of  source  reconstruction  for  real- 
world  activity  concentration  data  measured  by  an  operational 
network  of  sensors  (more  particularly  by  the  IMS  radionu¬ 
clide  network  that  is  maintained  under  the  auspices  of  the 
United  Nations  CTBT  Organization).  This  methodology  for 
source  reconstruction  has  been  applied  to  a  difficult  situation, 
namely,  reconstruction  of  the  source  characteristics  from  the 
CRL  release  using  only  a  small  number  of  activity  concen¬ 
tration  measurements  (8  measurements)  obtained  from  only 
three  sampling  stations.  The  problem  involved  utilization  of 
an  operational  backward  LS  model  for  long-range  transport 
by  the  atmosphere  on  a  continental  scale. 

The  principal  difficulty  in  the  reconstruction  lay  in 
providing  the  correct  a  priori  specification  of  the  model 
error  for  the  various  predicted  concentrations  associated  with 
the  measured  concentrations  used  for  the  source  parameter 
recovery.  A  naive  specification  for  the  model  error  gave  a 
source  reconstruction  that  was  quite  reasonable  in  terms  of 
the  accuracy  in  the  recovery  of  the  location  and  emission  rate, 
but  the  precision  in  the  estimates  was  generally  poor  in  the 
sense  that  the  reported  uncertainty  bounds  in  the  recovery  of 
the  source  parameters  did  not  include  the  actual  (true)  values 
for  these  parameters.  This  led  to  an  alternative  measurement 
model  which  incorporated  multipliers  (scale  factors)  with 
the  predicted  concentrations  in  order  to  compensate  for 
the  model  errors.  The  application  of  this  alternative  model 
was  shown  to  yield  significantly  improved  estimates  for  the 
source  parameters,  both  in  terms  of  their  accuracy  and  their 
precision. 

Although  Bayesian  probability  theory  offers  a  coherent 
and  rational  approach  for  source  reconstruction,  its  applica¬ 
tion  to  real-world  problems  using  real  sensor  networks  and 
operational  dispersion  models  will  require  a  better  under¬ 
standing  of  both  the  scale  and  structure  of  the  model  error  in 
the  predicted  concentrations.  It  is  anticipated  that  the  proper 
incorporation  of  information  about  the  underlying  model 
error  in  conjunction  with  Bayesian  inference  employing 


International  Scholarly  Research  Notices 


11 


m7 


- 1 - r 

- r- 

1 

1 

1 

1 

1 

-m— Tl  Lrm  iftir 

- 1 - 

r  irlTh-r^  ^ 

0  0.2  0.4  o.e 

m2 

i  0.8  1 

- 

1 1 
i 
i 
i 

i 

i 

n— n- n- TT"  dn-i-fTrSTT-nT-i-i-^  > 

0  0.2  0.4  0.6  0.8  1 

m4 

. . .  n  i-TITt  TTt|T 

1  r 

1 

1 

1 

1 

taim  i 

0  2  4 

6  8  10 

m6 

- 

i 

i 

i 

i 

1  i 

iuA 

0 - . n  hrnrfmoTTT^H 

0  2  4  6  8  10 

m8 


Figure  6:  One-dimensional  marginal  posterior  distribution  of  the  multipliers  m;  (/  =  1, 2, . . . ,  8)  used  to  compensate  for  the  model  error  in 
the  problem.  The  “actual”  values  of  the  multipliers  are  indicated  with  the  solid  vertical  line  which  should  be  compared  with  the  best  estimates 
of  the  multipliers  (posterior  means). 


state-of-the-science  MCMC  sampling  algorithms  (such  as  the 
MT-DREAM(ZS)  algorithm)  will  provide  the  most  flexible 
framework  for  source  reconstruction. 
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